Hands-on Exercise 2

First-order Spatial Point Patterns Analysis Methods

1.Overview

First-order Spatial Point Pattern Analysis (1st-SPPA) is a method used to study the spatial distribution and density of point-based phenomena—such as crime events, business locations, or public facilities—across a geographic area. Unlike higher-order analyses, 1st-SPPA focuses solely on individual point locations and their spread, without considering interactions between them. It helps identify where points are most concentrated, whether their distribution is uniform, and how widely they are dispersed. In this context, the chapter introduces hands-on techniques using the spatstat package to explore whether childcare centres in Singapore are randomly distributed, and if not, to pinpoint areas with higher concentrations.

2. Data

• Child Care Services data from data.gov.sg, a point feature data providing both location and attribute information of childcare centres.

• Master Plan 2019 Subzone Boundary (No Sea), a polygon feature data providing information of URA 2019 Master Plan Planning Subzone boundary data.

3. Installing and Loading R Packages

pacman::p_load(sf, terra, spatstat, 
               tmap, rvest, tidyverse)

4. Importing and Wrangling Geospatial Datasets

Import the master plan subzone boundary

mpsz_sf <- st_read("data/MasterPlan2019SubzoneBoundaryNoSeaKML (1).kml") %>% 
  st_zm(drop = TRUE, what = "ZM") %>% st_transform(crs = 3414)

Reading layer `URA_MP19_SUBZONE_NO_SEA_PL' from data source 
  `D:\ssinha8752\ISSS608-VAA\Hands-on_Ex\Hands-on_Ex02\data\MasterPlan2019SubzoneBoundaryNoSeaKML (1).kml' 
  using driver `KML'
Simple feature collection with 332 features and 2 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 103.6057 ymin: 1.158699 xmax: 104.0885 ymax: 1.470775
Geodetic CRS:  WGS 84

Extract the Key attributes and KML descriptions

extract_kml_field <- function(html_text, field_name) {
  if (is.na(html_text) || html_text == "") return(NA_character_)
  page <- read_html(html_text)
  rows <- page %>% html_elements("tr")
  value <- rows %>%
    keep(~ html_text2(html_element(.x, "th")) == field_name) %>%
    html_element("td") %>%
    html_text2()
  if (length(value) == 0) NA_character_ else value
}

Clean the structure and subzone data

mpsz_sf <- mpsz_sf %>%
  mutate(
    REGION_N = map_chr(Description, extract_kml_field, "REGION_N"),
    PLN_AREA_N = map_chr(Description, extract_kml_field, "PLN_AREA_N"),
    SUBZONE_N = map_chr(Description, extract_kml_field, "SUBZONE_N"),
    SUBZONE_C = map_chr(Description, extract_kml_field, "SUBZONE_C")
  ) %>%
  select(-Name, -Description) %>%
  relocate(geometry, .after = last_col())

Filter out Non-Mainland subzones

mpsz_cl <- mpsz_sf %>%
  filter(SUBZONE_N != "SOUTHERN GROUP",
         PLN_AREA_N != "WESTERN ISLANDS",
         PLN_AREA_N != "NORTH-EASTERN ISLANDS")

Save the cleaned subzone

write_rds(mpsz_cl, "data/mpsz_cl.rds")

Import the childcare services data

childcare_sf <- st_read("data/ChildCareServices.kml") %>% 
  st_zm(drop = TRUE, what = "ZM") %>%
  st_transform(crs = 3414)

Reading layer `CHILDCARE' from data source 
  `D:\ssinha8752\ISSS608-VAA\Hands-on_Ex\Hands-on_Ex02\data\ChildCareServices.kml' 
  using driver `KML'
Simple feature collection with 1925 features and 2 fields
Geometry type: POINT
Dimension:     XYZ
Bounding box:  xmin: 103.6878 ymin: 1.247759 xmax: 103.9897 ymax: 1.462134
z_range:       zmin: 0 zmax: 0
Geodetic CRS:  WGS 84

Validate the CRS Match

st_crs(mpsz_cl) == st_crs(childcare_sf)

[1] TRUE

4.1 Mapping Geospatial Datasets

DIY

library(sf)
library(ggplot2)
library(readr)

# Load cleaned subzone and childcare datasets
mpsz_cl <- read_rds("data/mpsz_cl.rds")
childcare_sf <- st_read("data/ChildCareServices.kml") %>%
  st_zm(drop = TRUE, what = "ZM") %>%
  st_transform(crs = st_crs(mpsz_cl))  # Ensure same CRS

Reading layer `CHILDCARE' from data source 
  `D:\ssinha8752\ISSS608-VAA\Hands-on_Ex\Hands-on_Ex02\data\ChildCareServices.kml' 
  using driver `KML'
Simple feature collection with 1925 features and 2 fields
Geometry type: POINT
Dimension:     XYZ
Bounding box:  xmin: 103.6878 ymin: 1.247759 xmax: 103.9897 ymax: 1.462134
z_range:       zmin: 0 zmax: 0
Geodetic CRS:  WGS 84

# Create the map
ggplot() +
  geom_sf(data = mpsz_cl, fill = "gray", color = "black", size = 0.2) +
  geom_sf(data = childcare_sf, color = "black", size = 1.2, alpha = 0.7) +
  labs(title = "Childcare Services Across Singapore Subzones",
       subtitle = "All layers aligned to EPSG:3414 (SVY21)",
       caption = "Source: URA Master Plan & Data.gov.sg") +
  theme_minimal()

Interactive mapping with tmap()

tmap_mode("view")

ℹ tmap mode set to "view".

tm_shape(childcare_sf) +
  tm_dots()

Registered S3 method overwritten by 'jsonify':
  method     from    
  print.json jsonlite

Switching back to plot

tmap_mode("plot")

ℹ tmap mode set to "plot".

5. Geospatial Data Wrangling

5.1 Converting sf data frames to ppp class

spatstat requires the point event data in ppp object form. The code chunk below uses [as.ppp()] of spatstat package to convert childcare_sf to ppp format.

childcare_ppp <- as.ppp(childcare_sf)

class(childcare_ppp)

[1] "ppp"

summary(childcare_ppp)

Marked planar point pattern:  1925 points
Average intensity 2.417323e-06 points per square unit

Coordinates are given to 11 decimal places

Mark variables: Name, Description
Summary:
     Name           Description       
 Length:1925        Length:1925       
 Class :character   Class :character  
 Mode  :character   Mode  :character  

Window: rectangle = [11810.03, 45404.24] x [25596.33, 49300.88] units
                    (33590 x 23700 units)
Window area = 796335000 square units

5.2 Creating owin object

it is a good practice to confine the analysis with a geographical area like Singapore boundary. In spatstat, an object called owin is specially designed to represent this polygonal region.

The code chunk below, as.owin() of spatstat is used to covert mpsz_sf into owin object of spatstat.

sg_owin <- as.owin(mpsz_cl)

class(sg_owin)

[1] "owin"

plot(sg_owin)